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Structural and dynamic properties of SPC/E water 

M. G. CampoP 

I have investigated the structural and dynamic properties of water by performing a series 
of molecular dynamic simulations in the range of temperatures from 213 K to 360 K, us- 
ing the Simple Point Charge- Extended (SPC/E) model. I performed isobaric-isothermal 
simulations (1 bar) of 1185 water molecules using the GROMACS package. I quantified 
the structural properties using the oxygen-oxygen radial distribution functions, order pa- 
rameters, and the hydrogen bond distribution functions, whereas, to analyze the dynamic 
properties I studied the behavior of the history-dependent bond correlation functions and 
the non-Gaussian parameter a.2(i) of the mean square displacement of water molecules. 
When the temperature decreases, the translational (r) and orientational (Q) order param- 
eters are linearly correlated, and both increase indicating an increasing structural order 
in the systems. The probability of occurrence of four hydrogen bonds and Q both have 
a reciprocal dependence with T, though the analysis of the hydrogen bond distributions 
permits to describe the changes in the dynamics and structure of water more reliably. 
Thus, an increase on the caging effect and the occurrence of long-time hydrogen bonds 
occur below ~ 293 K, in the range of temperatures in which predominates a four hydrogen 
bond structure in the system. 



I. Introduction 

Water is the subject of numerous studies due to its 
biological significance and its universal presence [1- 
3] . The thermodynamic behavior of water presents 
important differences compared with those of the 
other substances, and many of the characteristics of 
such behavior are often attributed to the existence 
of hydrogen bonds between water molecules. Scien- 
tists have found that the water structure produced 
by the hydrogen bonds is peculiar as compared to 
that of other liquids. Then, the advances in the 
knowledge of hydrogen bond behavior are crucial 
to understanding water properties. 
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The method of molecular dynamics (MD) allows 
to analyze the structure and dynamics of water at 
the microscopic level and hence to complement ex- 
perimental techniques in which these properties can 
be interpreted only in a qualitative way (infra-red 
absorption and Raman scattering [3], depolarized 
light scattering [5J E] , neutron scattering [7] , fem- 
tosecond spectroscopy [8-11] and other techniques 
[12-14]. 

Among the usual methods to study the short 
range order in MD simulations of water are the 
calculus of radial distribution functions, hydrogen 
bond distributions and order parameters. The ori- 
entational order parameter Q measures the ten- 
dency of the system to adopt a tetrahedral con- 
figuration considering the water oxygen atom as 
vertices of a tetrahedron, whereas the translational 
order parameter r quantifies the deviation of the 
pair correlation function from the uniform value of 
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unity seen in an ideal gas [HI \T§\ . The order pa- 
rameters are used to construct an order map, in 
which different states of a system are mapped onto 
a plane t-Q. The order parameters are, in general, 
independent, but they are linearly correlated in the 
region in which the water behaves anomalously |17j . 

The dynamics of water can by characterized by 
the bond lifetime, tub , associated to the process of 
rupturing and forming of hydrogen bonds between 
water molecules which occurs at very short time 
scale [9, 18, 19, 21-23]. t hb is obtained in MD 
using the history-dependent bond correlation func- 
tion P(i), which represents the probability that an 
hydrogen bond formed at time t = remained con- 
tinuously unbroken and breaks at time t [241 125] . 

Also, the dynamics of water can be studied by 
analyzing the mean-square displacement time se- 
ries M{t). In addition to the diffusion coefficient 
calculation at long times in which M(t) oc t, in the 
supercooled region of temperatures and at interme- 
diate times M(i) oc t a (0 < a < 1). This behavior 
of M(i) is associated to the subdiffusive movement 
of the water molecules, caused by the caging effect 
in which a water molecule is temporarily trapped 
by its neighbors and then moves in short bursts due 
to nearby cooperative motion. A time t* character- 
izes this caging effect (see Sec. II for more details) 

In a previous work, we found a g-exponential be- 
havior in -P(i), in which q increases with T -1 ap- 
proximately below 300 K. q(T) is also correlated 
with the probability of occurrence of four hydrogen 
bonds, and the subdiffusive motion of the water 
molecules [25] . 

The relationship between dynamics and struc- 
tural properties of water has not been clearly es- 
tablished to date. In this paper, I explore whether 
the effect that temperature has on the water dy- 
namics reflects a more general connection between 
the structure and the dynamics of this substance. 

II. Theory and method 

I have performed molecular dynamic simulations of 
SPC/E water model using the GROMACS package 
[13 [3D], simulating fourteen similar systems of 1185 
molecules at 1 bar of pressure in a range of temper- 
atures from 213 K to 360 K. I initialized the system 
at 360 K using an aleatory configuration of water 



Table 1: Details of the simulation procedure. Du- 
ration of the stabilization period (t es t) and the MD 
sampling (Jm d ) in the different ranges of tempera- 
tures 



Temp, range (K) 


test (ns) 


t M D (ns) 


213 - 243 


20.0 


10.0 


253 - 273 


16.0 


10.0 


283 - 360 


16.0 


8.0 



molecules, assigning velocities to the molecules ac- 
cording to a Boltzmann's distribution at this tem- 
perature. For stabilization, I applied Berendsen's 
thermal and hydrostatic baths at the same tem- 
perature and 1 bar of pressure [31 j - Then, I ran 
an additional MD obtaining an isobaric-isothermal 
ensemble. I obtained the other systems in a similar 
procedure, but using as initial configuration that 
of the system of the preceding higher temperature 
and cooling it at the slow rate of 30 K ns -1 [17] . 
Stabilization and sampling periods for the systems 
at different temperatures are indicated in Table [TJ 
Simulation and sampling time steps were 2 fs and 
10 fs, respectively. The sampling time step was 
shorter than the typical time during which a hydro- 
gen bond can be destroyed by libration movements. 

I calculated the hydrogen bond distribution func- 
tions f(n) (n = 0,1,..., 5), which is the probabil- 
ity of occurrence of n hydrogen bonds by molecule, 
considering a geometric definition of hydrogen bond 
[20] . As parameters for this calculation, I used a 
maximum distance between oxygen atoms of 3.5 A 
and a minimum angle between the atoms Odonor~ 
H-O accep t or of 145°. 

The radial distribution function (RDF) is a stan- 
dard tool used in experiments, theories, and sim- 
ulations to characterize the structure of condensed 
matter. Using RDFs, I obtained the average num- 
ber, N, of water molecules in the first hydration 
layer (the hydration number) 

N = Airp f g(r)r 2 dr (1) 
Jo 

where p is the number density. 
The translational order parameter, r, is defined 
in Ref. [TB] as 
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Figure 1: Hydrogen bonds distribution functions 
f(n) (n = 0, 5) versus T. The zones A, B and C 
correspond to ranges of temperatures in which oc- 
cur different relationships between /(4),/(3) and 
/(2). Note the reciprocal scale for the tempera- 
tures. See the text for details. 
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Figure 2: Oxygen-oxygen radial distribution func- 
tions for the systems at 213 K (continuous line), 
293 K (dashed line), and 360 K (dotted line). In- 
set: The hydration number N vs. T 5 . 



\g(s)-l\ds 



(2) 



where the dimensionless variable s = rn 1 / 3 is the 
radial distance r scaled by the mean intermolecular 



distance n 1 / 3 , and S c corresponds to half of the 
simulation box size. 

The orientational order parameter Q is defined 
as UM 



I o N 4 4 
\ i=l j=l l=j+l 



cosQ 



jik 



(3) 



where 6 jik is the angle formed by the atoms Oj— 0$- 
Ofc. Here, is the reference oxygen atom, and Oj 
and Ok are two of its four nearest neighbors. Q=l 
in an ideal configuration in which the oxygen atoms 
would be located in the vertices of a tetrahedron. 

I obtained the bond correlation function P(t) 
from the simulations by building a histogram of 
the hydrogen bonds lifetimes for each configura- 
tion. Then, I fitted this function with a Tsallis 
distribution of the form 



exp, 



(4) 



being t the hydrogen bond lifetime and q the nonex- 
tensivity parameter [551 132] ■ If t/ = 1, Eq. @D 
reduces to an exponential, whereas if q > 1, P(t) 
decays more slowly than an exponential. This last 
behavior occurs when long lasting hydrogen bonds 
increase their frequency of occurrence. 

The subdiffusive movement of water occurs when 
the displacement of the molecules obeys a non- 
Gaussian statistics. This behavior is characterized 
by t* , the time in which the non-Gaussian param- 
eter at2(t) reaches a maximum [see Eq. (J5J) ] . Then, 
t* is the parameter associated to the average time 
during which a water molecule is trapped by its en- 
vironment (caging effect), and this prevents it from 
reaching the diffusive state [57] . 



a 2 (t) 



1 



(5) 



III. Results and discussion 



Three zones or ranges of temperatures can be dis- 
tinguished in the graph of the hydrogen bond dis- 
tributions f(n) vs. T (see Fig. [TJ. Zone A (T > 
350 K) in which /(3) > /(2) > /(4), zone B (293 
K > T > 350 K) in which /(3) > /(4) > /(2), and 
zone C (T < 293 K) in which /(4) > /(3) > /(2). 
These results indicate a predominant structure of 
three and two hydrogen bonds (3HB-2HB) in zone 
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A, 3HB-4HB in zone B, and 4HB-3HB in zone C, 
respectively. /(4) oc T^ 1 in all ranges of tempera- 
tures, showing that the tetrahedral structure of wa- 
ter decreases with the increase of this variable. /(3) 
increases with T up to 293 K, and then remains 
approximately constant (~ 0.4) up to 360 K. f(2) 
also increases with T in all range of temperatures, 
but only overcomes /(4) at T > 350 K. 
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Figure 3: Position of the first minimum of the 
oxygen-oxygen radial distribution function vs T 4 , 
associated to the size of the first hydration layer. 

Fig. [5] shows the oxygen-oxygen RDFs corre- 
sponding to the systems at 213 K, 293 K and 360 
K. When the temperature decreases, the minimum 
and maximum tend to be more defined. This being 
associated with an increasing order in the system. 
The position of the first minimum moves closer to 
the origin decreasing the size of the first hydration 
layer (oc T A see Fig. [3]). Both facts can be associ- 
ated to the decrease of the hydration number from 
N ~ 5 to N ~ 4 (see inset, Fig. [2J. 

The simultaneous behavior of Q and r is shown 
in the order map of Fig. HJ in which the location of 
the values corresponding to 293 K are indicated by 
an arrow. The order parameters present similar be- 
haviors with the temperature. Upon cooling, these 
parameters are linearly correlated and move in the 
order map along a line of increasing values, up to 
reaching maximum values at 213 K. The slope of 
the line increases a little for T > 293 K, indicating 
that r has a response to the increase of T slightly 
higher than Q in this range of temperatures. The 




Figure 4: Order map with the values of the or- 
der parameters corresponding to the simulated sys- 
tems. Note the change in the slope of the line at 
T ~ 273 K 



positive values of the slopes indicate an increasing 
order of the system when the temperature decrease. 

The f(ri) functions allow to obtain a more 
detailed picture of the structural orientational 
changes at shorter ranges between water molecules 
than the orientational order parameter. While a 
small change at 293 K occurs in the order map, the 
structures of two, three and four hydrogen bonds 
are alternated in importance when the temperature 
changes. The ability of f(n) to more reliably de- 
scribe the structure of the water occurs because the 
calculation of the hydrogen bond distributions in- 
cludes the location of the hydrogen atoms, whereas 
Q only quantifies the changes in the average an- 
gle between neighbor oxygen atoms. Although the 
behavior of /(4) and Q are correlated (see Fig. 
[SJ, /(4) shows a greater response to the temper- 
ature than Q, indicating that the main change in 
the tetrahedral structure with the decrease of the 
temperature occurs mainly in the orientation of 
the bonds between water molecules. The approx- 
imately linear correlation between both variables 
also indicates a similar dependence with the tem- 
perature (oc T -1 ). 

Figure [5] shows the behavior of the dynamical pa- 
rameters t* and q with Q. The characteristic time 
t* has an exponential response to Q > 0.58 (T < 
360 K), but the slope of the semilog plot of t* vs. 
Q increases significantly for Q > 0.67. A similar 
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Figure 5: Q vs /(4). The change in /(4) is higher 
than that of Q in the range of temperatures stud- 
ied. See the text for details. 

change occurs for Q > 0.67 in the linear correla- 
tion between q and Q. Then, the values Q ~ 0.67 
and r 1.1 of the order map can be associated 
to changes in the dynamics of the system. The 
transition of q ss 1 to q > 1 indicates the increase 
of the probability of two water molecules remain- 
ing bonded by a hydrogen bond during an unusual 
long time, whereas the increase of t* is associated to 
the increase of the time during which the molecules 
remain in a subdiffusive regime. 

However, only the analysis of the f(n) functions 
reveals the structural modification that explains 
the structural and dynamic changes in the system. 
The changes in the increase of the order map, t* (Q) 
and q(Q) occur below 293 K, in the range of tem- 
peratures in which prevail a structure of four hy- 
drogen bonds in the system. 

IV. Conclusions 
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Figure 6: (a) Semilog plot of t* vs. Q. (b) q vs. Q. 
Sec the text for details. 



of occurrence of four hydrogen bonds and the ori- 
entational order parameter Q. However, only the 
analysis of the behavior of the hydrogen bond dis- 
tribution allows to deduce that, when a tetrahedral 
structure associated to the percentage of four hy- 
drogen bonds predominates, the behavior of the dy- 
namical variables P{t) and t* show the occurrence 
of long lasting hydrogen bonds and caging effect 
between the molecules of the system. 



The molecular dynamic method allows to study the 
structure and dynamics of the SPC/E model of wa- 
ter in the range of temperatures from 213 K to 360 
K. 

Lowering the temperature of the system from 
360 K to 213 K, the number of water molecules 
in the first hydration layer decreases from TV <~ 5 
to TV ~ 4, along with a decrease in size. The in- 
crease of the tetrahedral structure of the system is 
also characterized by a growth of the percentage 
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